Density Functionals that Recognize Covalent, Metallic, and Weak Bonds 



Jianwei Sun, 1 Bing Xiao, 1 Yuan Fang, 1 Robin Haunschild, 2 Pan Hao, 1 Adrienn 
Ruzsinszky, 1 Gabor I. Csonka, 3 Gustavo E. Scuseria, 2 ' ' 5 and John P. Perdew 1 

1 Department of Physics and Engineering Physics and Quantum Theory Group, 
Tulane University, New Orleans, Louisiana 70118, USA 
2 Department of Chemistry, Rice University, Houston, Texas 77005, USA 
3 Department of Inorganic and Analytical Chemistry, 
Budapest University of Technology and Economics, H-1521 Budapest, HUNGARY 
J ' Department of Physics and Astronomy, Rice University, Houston, Texas 77005, USA 
Chemistry Department, Faculty of Science, King Abdulaziz University, Jeddah 21589, Saudi Arabia 

(Dated: June 4, 2013) 

Computationally-efficient semilocal approximations of density functional theory at the level of 
the local spin density approximation (LSDA) or generalized gradient approximation (GGA) poorly 
describe weak interactions. We show improved descriptions for weak bonds (without loss of ac- 
curacy for strong ones) from a newly-developed semilocal meta-GGA (MGGA), by applying it to 
molecules, surfaces, and solids. We argue that this improvement comes from using the right MGGA 
dimensionless ingredient to recognize all types of orbital overlap. 
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Due to its computational efficiency and reasonable 
accuracy, the Kohn-Sham density functional theory 
[IH2] with semilocal approximations to the exchange- 
correlation energy, e.g., the local spin density approxi- 
mation (LSDA) [H [5] and the standard Perdew-Burkc- 
Ernzcrhof (PBE) generalized gradient approximation 
(GGA) [6 ], is one of the most widely- used electronic 
structure methods in materials science, surface science, 
condensed matter physics, and chemistry. Semilocal ap- 
proximations display a well-understood error cancella- 
tion between exchange and correlation in bonding re- 
gions. Thus some intermediate-range correlation effects, 
important for strong and weak bonds, are carried by 
the exchange part of the approximation. However, it is 
well-known that these approximations cannot yield cor- 
rect long-range asymptotic dispersion forces [7j. This 
raises doubts about the suitability of semilocal approxi- 
mations for the description of weak interactions (includ- 
ing hydrogen bonds and van der Waals interactions), even 
near equilibrium where most interesting properties occur. 
These doubts are supported by the performance of LSDA 
and GGAs, which are not very useful for many important 
systems and properties (such as DNA, physisorption on 
surfaces, most biochemistry, etc.). 

However, these doubts are challenged by recent devel- 
opments in semilocal meta-GGAs (MGGA) [SHU] (which 
are useful by themselves and as ingredients of hybrid 
functionals [33]). Compared to GGAs, which use the 
density n(r) and its gradient Vn as inputs, MGGAs 
additionally include the positive kinetic energy density 
t = J2k I^V'fcl 2 /2 of the occupied orbitals Vfc- F° r 
simplicity, we suppress the spin here. By including 
training sets of noncovalent interactions, the molecule- 
oriented and heavily-parameterized M06L MGGA was 
trained to capture medium-range exchange and corre- 
lation energies that dominate equilibrium structures of 
noncovalent complexes [SJ. Madsen et al. showed that 



the inclusion of the kinetic energy densities enables MG- 
GAs to discriminate between dispersive and covalent in- 
teractions, which makes the M06L MGGA [9J suitable 
for layered materials bonded by van der Waals inter- 
actions [TSJ [IB]. Besides improvement for noncovalent 
bonds, simultaneous improvement for metallic and cova- 
lent bonds is also an outstanding problem for semilocal 
functionals [171 HE]- Pef. [T5] has shown that the revised 
Tao-Perdew-Staroverov-Scuseria (revTPSS) [TU] MGGA, 
due to the inclusion of the kinetic energy density, simul- 
taneously predicts accurate results for the adsorption en- 
ergy of CO on the Pt (111) surface and the lattice con- 
stant and surface energy of the substrate, while GGA and 
LSDA jl8] do not. These successful applications show 
the advantages and flexibility brought by inputting the 
kinetic energy density for MGGAs. The effect of the 
dependence on the kinetic energy density has also been 
studied recently in Refs.fTT | [T4 t leading to a new MGGA 
called MGGA_MS2 |14| . However, an important question 
for the development of MGGAs and thus for understand- 
ing the results from MGGAs remains unanswered: how 
should the kinetic energy density be built into MGGAs, 
why, and what are the consequences? 

To answer this question, let's look at the input, the 
positive kinetic energy density, and how it is used in the 
M06L, revTPSS, and MGGA.MS2 MGGAs. r can be 
expressed explicitly in terms of the density n(r) for two 
kinds of paradigm systems. For one- and two-electron 
ground states, r reduces to the von Weizsacker kinetic 
energy density, t w — ||Vn| 2 /n. These systems, espe- 
cially the hydrogen atom, are paradigmatic for quantum 
chemistry. On the other hand, for a uniform electron gas 
of density n, which is a paradigm density for solids and 
a constraining limit for most approximate functionals, r 
becomes r unlf = (37r 2 ) 2 / 3 n 5 / 3 . In the uniform electron 
gas, electrons are fully delocalized and orbitals are highly 
overlapped. 
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TABLE I: Values of z, a, and t 1 for typical regions. 



Region 


z a 


t- 1 


Single orbital 


1 


5s 2 /3 


Slowly- varying density 




« 1 


Overlap of closed shells 


w0 > 1 


> 1 



With t w and r unlf , three different dimensionless pa- 
rameters can be constructed: z = t w /t, a — (r — 
T w )/r unil , and t" 1 = r/r unif = a + 5s 2 /3. Here, 
s = |Vn|/[2(37r 2 ) 1 / 3 n 4 / 3 ] is the reduced density gradient, 
an important dimensionless parameter measuring the in- 
homogeneity of the density and employed in GGAs. Ta- 
ble [I] shows the values of these three parameters for three 
different typical regions: 1) regions of one- and two- 
electron densities that characterize single bonds; 2) re- 
gions of slowly- varying density that characterize metallic 
bonds; and 3) regions of density overlap between closed 
shells that characterize noncovalent bonds. It can be 
seen that only a can recognize all types of orbital overlap 
and should be chosen as the ingredient for constructing 
MGGAs. It equals for the regions of one- and two- 
electron densities and « 1 for regions of slowly-varying 
density. For regions of density overlap between closed 
shells, t w w because the density gradient is small by 
symmetry, r/r unlf oc n/n 5 / 3 , and thus a 3> 1 as n is 
small in such regions. Note a > since t w is a lower 
bound on r [19]. 

However, instead of a, which is only used for recover- 
ing the fourth-order gradient expansion of the exchange 
energy of a slowly- varying density, revTPSS uses z to 
identify different orbital-overlap regions, z identifies the 
single-orbital region with z = 1 and the slowly-varying 
density with z « 0, but it cannot distinguish regions of 
slowly- varying density from those of density overlap be- 
tween closed shells, both of which have z ~ 0. Since 
revTPSS used the information of the hydrogen atom 
(where z = 1) and the slowly- varying density (where 
z ~ 0) to fix its parameters and orient toward cova- 
lent bonds, this limitation of z makes revTPSS unable 
to track noncovalent bonds. This is manifested by its 
performance on the problem of graphene adsorbed on 
the Ni(lll) surface. 

Graphene was found to bind to the Ni (111) surface 
at 2.11±0.07 A chemically (from experiment [20]) and 
at 3.3 A physically (from theory [3JJ), respectively. In 
Fig. l(a)[ the binding curve of the random phase approx- 
imation (RPA) from Ref [HH which naturally includes 
the nonlocal van der Waals interaction, is included as 
the reference. RPA predicts both the chemisorption and 
the physisorption minima with the correct long-range 
asymptotic behavior [21]. Note that RPA calculations 
are about 2 orders of magnitude more time consuming 
than those of semilocal functionals. 



the physisorption one, while PBE misses the latter and 
almost misses the former. Compared to the RPA results, 
revTPSS predicts the chemisorption with a better bind- 
ing energy than those of LSDA and PBE. But it still 
misses the physisorption minimum, consistent with the 
result of Ref. [22] indicating its inability to capture nonco- 
valent bonds. Although all semilocal functionals consid- 
ered here yield wrong long-range asymptotic behaviors, 
MGGAJV1S2 and M06L remarkably capture the double 
minima at accurate distances, demonstrating their ability 
to describe both strong and weak bonds. If we use the ex- 
perimental and RPA binding distances (2.11±0.07 A, 3.3 
A) as references for the chemisorption and physisorption, 
respectively, then MGGA_MS2 (2.09 A, 3.29 A) yields 
a better agreement than M06L (2.29 A, 3.25 A) [To] . 
In terms of the binding energies, M06L (64 meV, 64 
meV) [16, is better than MGGAJVIS2 (49 meV, 23 meV) 
when compared to RPA (67 meV, 60 meV) [21] , with the 
first number in the parentheses for the chemisorption and 
the second for the physisorption. 

In MGGAJV1S2 [14]. a is the only ingredient built from 
r, and it can recognize all types of orbital overlap. Al- 
though the two fitting parameters of MGGA_MS2 are 
determined by two data sets that do not involve weak 
interactions [14], the ability of a to identify the over- 
lap between closed shells helps MGGA_MS2 to capture 
weak interactions near equilibrium, which explains why 
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shows that 

LSDA predicts the chemisorption minimum but misses 
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MGGAJV1S2 captures both minima in Fig. 1(a) 

M06L, on the other hand, uses i _1 to identify different 
orbital-overlap, with a only to make the one-electron cor- 
relation energy exact. Table [i] shows that t^ 1 ss 1 for the 
slowly- varying density, and t~ x 3> 1 for the overlap be- 
tween two closed shells. Madsen et al. have shown that 
i _1 can discriminate between covalent and noncovalent 
interactions [T5] , which then explains the ability of M06L 
to capture the double minima in Fig. 1(a) This exam- 
ple together with many other successful applications of 
M06L, including CO adsorbed on the Pt(lll) surface [H] 
and layered solids [TS], seems to suggest that t _1 , in ad- 
dition to s, would be a good dimensionless parameter 
that measures the inhomogeneity of the density. How- 
ever, t~ 
regions 

tion in comparison with s. It can be seen in Fig. 1(b) 
that, for the region around the nucleus of Ar where the 
Is orbital dominates, t _1 resembles s in shape. To over- 
come the shortcoming of t^ 1 being unable to identify 
single-orbital regions, and to fit to a large number of 
data sets that include covalent single bonds (a » 0), 
35 fitting parameters had to be introduced into M06L, 
causing an oscillation in the exchange enhancement fac- 
tor near small a as shown in Fig. 1(c) The exchange 
enhancement factor F x characterizes the enhancement of 
the exchange energy density with respect to its local ap- 
proximation, defined by E x [n\ — J d 3 rne x nli (n)F x (s, a), 
where e" nlf (n) = — 3(37r 2 n) 1 / 3 /47r is the exchange energy 
per electron of the uniform electron gas. 

The oscillation in the M06L F x can cause two major 



1 loses the ability to identify the single-orbital 
where = 5s 2 /3 provides no extra informa- 
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FIG. 1: (a) Binding-energy curves for graphene adsorption on the Ni (111) surface. The M06L and RPA values are from 
Refs. fT6l and 1211 respectively, b) Distributions of s, a, z, and i _1 for Ai2 at equilibrium, (c) Exchange enhancement factor, 
(d) Binding-energy curve of Ar2 for different functionals. See Ref. 1231 for computational details. 



problems. The first one is the lack of computational sta- 
bility and the high requirements to converge calculations. 
An example is the oscillation of the binding curve of the 
Ar dimer shown in Ref 23, The oscillation can be re- 
moved by using much higher computational settings for 
the calculations, but even then M06L still yields a too- 
shallow minimum far away from the experimental bind- 
ing distance, as shown in Fig. 1(d) The performance 



of M06L on the Ar dimer and its inability to bind the 
Kr 2 and Xe 2 rare gas dimers |25j . which are typical sim- 
ple tests for van der Waals interactions, contradict its 
good performance on many other noncovalent systems. 
This also leads to the second problem, the consistency in 
describing different bonds in different chemical environ- 
ments. To investigate this, we study the performance of 
M06L on systems that are far outside its training sets, 
e.g., lattice constants of solids. We choose the six ionic 
insulators and semiconductors studied in Ref [5^1 Ref l2"rJI 
showed that, although the covalent and ionic bonds domi- 
nate, the van der Waals interaction plays an essential role 
in determining accurate lattice constants. It found that 
the errors in lattice constants of these six solids are re- 
duced by a factor of 2, in comparison with experimental 
data, when the van der Waals interactions were included 
on top of the PBE functional. 

Tabic |n] shows that LDA undcrstimates the lattice 
constants of the 6 solids while PBE overestimates them. 
revTPSS improves the lattice constants over PBE slightly 



for ionic insulators but significantly for semiconduc- 
tors when compared to experimental data. However, 
revTPSS still yields too large lattice constants because 
revTPSS misses weak interactions in these solids. There- 
fore, the improvement of revTPSS over PBE should 
be ascribed to its better description on covalent sin- 
gle bonds, by which those semiconductors are bound. 
Adding the van der Waals interactions on top of revTPSS 
would be expected to further improve the lattice con- 
stants. This is largely realized by MGGA_MS2 alone, 
which reduces the mean absolute error (MAE) from 0.042 
A of revTPSS to a remarkable 0.011 A. This demon- 
strates one more time that MGGA_MS2 describes well 
both strong (covalent and ionic) bonds and weak interac- 
tions because a can recognize all types of orbital overlap. 

The pattern of the M06L performance on these 6 solids 
however is inconsistent. Compared to experimental data, 
M06L predicts accurate lattice constants for C and Si. 
But it yields too-large ones, even larger than those of 
PBE, for Ge and GaAs which have relatively large ionic 
radii and thus significant van der Waals interactions. 
Similar behavior is also observed for the two ionic solids. 
M06L yields an accurate lattice constant for MgO, but 
a too-large one for NaCl. This random performance 
is likely related to the oscillation of the exchange en- 
hancement factor near small a and small s as shown in 



Fig. 1(c) because small s is more relevant to solids than 



to molecules. The inconsistent performance of M06L in 
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this test set cautions against its use in solids, and is 
rooted in using t^ 1 to construct the functional. 



TABLE II: Errors in lattice constants of ionic insulators and 
semiconductors from different functionals. The zero-point an- 
harmonic expansion has been removed from the experimental 
lattice constants [271. ME: mean error; MAE: mean absolute 
error. 



solids 


LDA PBE 


M06L 


revTPSS MGGA_MS2 Expt. 


C 


-0.022 0.0f4 


-0.010 


0.003 


-0.007 


3.555 


Si 


-0.0f7 0.046 


-0.010 


0.017 


0.005 


5.422 


Ge 


-0.013 0.124 


0.137 


0.038 


0.005 


5.644 


GaAs 


-0.026 0.111 


0.143 


0.039 


0.002 


5.641 


NaCl 


-0.098 0.130 


0.117 


0.102 


0.029 


5.565 


MgO 


-0.018 0.073 


0.012 


0.052 


0.019 


4.188 


ME 


-0.032 0.083 


0.065 


0.042 


0.009 




MAE 


0.032 0.083 


0.072 


0.042 


0.011 





The above discussions qualitatively explain the perfor- 
mance of revTPSS, M06L, and MGGA_MS2, with a spe- 
cial focus on noncovalent interactions, whereas quantita- 
tive analysis using enhancement factors is a subtle issue. 
Madsen et al. [15] argued that the inability of TPSS [8] 
(whence also of revTPSS) to capture weak interactions is 
due to the increase of its F x at small s as a becomes large, 
as shown in Fig. 1(c) for revTPSS. That argument is sup- 
ported by the fact that the MGGAJMS2 F x at small s 



decreases monotonically as a becomes large; see Fig. |l(c) 
of this work and Fig. 2 of Ref. 11} The monotonically 
decreasing a-dependence can be rationalized by the fact 
that the strength of different chemical bonds decreases 
with increasing a (e.g., a — for single bonds, a ~ 1 for 
metallic bonds, and a 1 for noncovalent bonds.) Note 
that a is related to the electron localization function [28] , 
r/=l/(l+a 2 ), which has been used to establish a rigor- 
ous topological classification of chemical bonds [2S], but 
surprisingly was not recognized as a useful ingredient for 
MGGAs. However, it should be stressed that the im- 
provement for weak interactions near equilibrium from 
MGGA_MS2 results from the good balance between its 
s— and a— dependences [i"4"l |2"3"] . 

To show its potential impact on a wide range of re- 
search and also to test its consistency for different sys- 
tems, we apply MGGA_MS2 for stacking energies of nu- 
cleobases, which are essential for the structural stability 
and function of DNA and RNA (30] [31]. Fig.[2]shows how 
the interaction energies for several stacked pairs of nu- 
cleobases (from adenine (A), cytosine (C), guanine (G), 
and uracil (U)) vary as one nucleobase rotates relative to 
the other around their perpendicular axis. The reference 
point for the rotation was defined using the method of 
Elstner et al. [35] , in which the center of mass of each 
nucleobase was aligned with the glycocidic bonds paral- 
lel. A counterclockwise (right-handed) rotation was then 
applied (see Fig. 2 of Ref. 31). In each case, the nu- 
cleobase structure was fixed to that of the optimized 
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FIG. 2: Interaction energy (in kcal/mol) as a function of ro- 
tation angle (in degrees) for stacked homonucleobase dimers. 
The MP2 values were calculated at a fixed separation distance 
of 3.4 A, taken from Ref. [31] The PBE, revTPSS, M06L, and 
MGGA_MS2 values were calculated at their own optimized 
distances for 180° . See Ref. [53] for computational details. 1 
kcal/mol = 0.0434 eV. 



isolated molecule. Except for the CC stacking struc- 
tures, where MGGA_MS2 misses the minimum at 180° 
of MP2 (M0ller-Plesset second-order perturbation theory 
in the electron-electron interaction), MGGAJMS2 yields 
the stacking energies qualitatively in agreement with the 
reference MP2 values [3T] . For all considered cases where 
binding between nucleobases is expected, MGGAJMS2 
significantly improves the stacking energies over PBE 
and revTPSS, but still underestimates them about the 
same amount as M06L overestimates, in comparison 
with the reference MP2 method. If this performance 
persists for other pairs, then MGGA_MS2 can provide 
a better description for DNA and RNA conformations 
than PBE (and, when long-range van der Waals interac- 
tions are added [35], probably better than M06L). This 
could be an important step toward computer-assisted 
drug design since the computational cost of the semilocal 
MGGA_MS2 is affordable for large bio-molecules. 

The computationally-efficient semilocal density func- 
tionals must fail for many stretched bonds |55j . But we 
have suggested that some MGGAs can be usefully accu- 
rate for unstretched or modestly-stretched ones, not just 
for strong but even for weak bonds. Successful but simple 
MGGAs employ only a = (r — T W )/r umf to recognize all 
types of orbital overlap, and extrapolate monotonically 
from < a < 1 to large a. This insight should guide 
the construction of further-improved semi-empirical and 
nonempirical density functionals. 
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